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Elucidating the fitness measures optimized during the evolution of complex biological systems 
is a major challenge in evolutionary theory. We present experimental evidence and an analytical 
framework demonstrating how biochemical networks exploit optimal control strategies in their evo- 
lutionary dynamics. Optimal control theory explains a striking pattern of extremization in the 
redox potentials of electron transport proteins, assuming only that their fitness measure is a control 
objective functional with bounded controls. 
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I. INTRODUCTION 

In a famous paper sent to Charles Darwin and pre- 
sented before the Linnean Society in 1858 [l], Alfred Rus- 
sel Wallace - often considered the co-discoverer of natural 
selection - proposed that evolution exploits the principles 
of feedback control in the generation of biological com- 
plexity. Wallace stated: 



"The action of this [evolutionary selection] principle is 
exactly like that of the centrifugal governor of the steam 
engine, which checks and corrects any irregularities al- 
most before they become evident; and in like manner no 
unbalanced deficiency in the animal kingdom can ever 
reach any conspicuous magnitude, because it would make 
itself felt at the very first step, by rendering existence 
difficult and extinction almost sure soon to follow" [3. 



During the ensuing development of evolutionary the- 
ory, the possibility that nature employs evolutionary con- 
trol strategies to maximize the fitness of biological net- 
works has often been discussed in the context of cybernet- 
ics, the study of self-regulation. However, to our knowl- 
edge, no direct, quantitative evidence for Wallace's con- 
tention - namely, that evolutionary dynamics itself may 
be self-regulating - has ever been reported. In this pa- 
per, we provide such evidence, and develop a quantitative 
physical theory for the interrogation of control phenom- 
ena in the evolution of biological systems. 

Evolution is guided by the optimization of fitness mea- 
sures that balance functionally beneficial properties. In 
modern theories of evolutionary dynamics, such as the 
quasispecies model Q and variants thereof, the fitness 
measure of a biological system plays a role analogous 
to that of the free energy of a mechanical system. The 
dynamics of the system, embodied through mutations, 
seeks to optimize this measure. Recently, with advances 
in the understanding of molecular biophysics, increasing 
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attention has been paid to characterizing the fitness mea- 
sures underlying the evolution of proteins. For example, 
simulations of protein sequence evolution have confirmed 
that protein cores evolve almost universally to maximize 
the free energy gap between the folded and denatured 
states 0. However, for functional properties of proteins 
and protein networks, the appropriate biological fitness 
measures are not so clear Q • A current challenge in evo- 
lutionary theory is to identify how the fitness measures of 
complex biological systems depend on the physical prop- 
erties of their constituent proteins. 

In the hierarchical evolution of protein networks, bio- 
logical self-organization Q influences the dynamics that 
occur on shorter time scales. Although most theories of 
evolutionary dynamics have modeled evolution as a dy- 
namical system seeking to optimize a potential or free en- 
ergy, multi-timescale evolution of protein networks may 
be modeled within a broader framework as a control 
problem. Optimal control (OC) theory is generally con- 
cerned with the determination of the time-dependent 
functional form of the Hamiltonian of a controlled dy- 
namical system that maximizes a desired objective func- 
tion j6|. An important difference between a dynamical 
system and a control system is that the latter distin- 
guishes between the free dynamics of the system and 
the dynamics regulated by controls. In the present case, 
these controls can take the form of functional protein 
properties. 

The evolution of a biological system may be modeled as 
a control system if the regulatory functional properties 
of its constituent proteins coevolve with the network's 
overall function. Should the evolutionary dynamics of 
such a system demonstrate features indicative of opti- 
mal controls, this would constitute evidence that the 
system's evolution has attained a sophisticated level of 
self-organization amounting to the solution of an OC 
problem. Here, we show that application of this theory 
to active site mutations in an enzyme network of central 
importance for metabolism - the electron transport chain 
(ETC) 7] - indicates that the redox potentials of electron 
transport proteins are controlling the evolutionary dy- 
namics of this network in an optimal fashion, providing 
insight into the self-organization of this system. 
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II. ARTIFICIAL EVOLUTION OF ELECTRON 
TRANSPORT PROTEINS 

The mitochondrial electron transport chain removes 
electrons from the high-energy electron donor NADH and 
passes them to the electron acceptor O2 through a series 
of redox reactions involving electron transport proteins. 
These reactions are coupled to the generation of a pro- 
ton concentration gradient across the mitochondrial in- 
ner membrane, which is ultimately used to produce ATP. 

Several of these electron transport proteins (for ex- 
ample, NADH-Q reductase, cytochrome reductase and 
cytochrome c oxidase) act as both electron carriers and 
proton pumps, simultaneously catalyzing the transfer of 
protons against their concentration gradient. Although 
the molecular mechanisms of these proton pumps were 
unclear for some time 0, [H, recent work has begun to 
explore these mechanisms, particularly that of the termi- 
nal protein in the chain, cytochrome c oxidase [lol . [llj | . 
Belevich et al. [Ill] studied the proton pump mechanism 
of this protein in real-time by spectroscopic and electro- 
metric techniques after laser-activated electron injection 
into the oxidized enzyme. It was found that the electron 
transfer reaction to the primary heme site ("heme a") 
of this protein raises the pKa of a proton "pump site" 
amino acid side chain due to the proximity of the associ- 
ated negative charge; importantly, the pump protonation 
site is not in the immediate vicinity of heme a, and its 
distance to the heme a plays a role in determining the 
thermodynamic efficiency of energy transduction. The 
increased pKa of the pump site draws a proton from the 
interior of the mitochondrion, with the catalytic assis- 
tance of several amino acid residues including Aspl24, 
Glu278, and Lys354 in P. denitrificans Protona- 
tion of the pump site initiates additional protonation and 
electron transfer reactions involving distinct sites distal 
to heme a. In particular, a second protonation induces 
release of the "pump" proton outside the mitchondrial 
membrane due to electrostatic repulsion, thus increasing 
the proton concentration gradient. Note that alteration 
of the redox potential of heme a would alter not only 
the thermodynamic efficiency of energy transduction but 
also the kinetics of its catalysis by Aspl24, Glu278 and 
Lys354, such that mutations around these sites distant 
from heme a would be required to maximize the efficiency 
of energy transduction. In extreme cases, mutations that 
alter the redox potential may even render energy trans- 
duction impossible J^]. 

Our prior work (isl . IT^ explored the mapping between 
redox potential and amino acid sequence in the heme 
microenvironment of ETC proteins. We pursued a strat- 
egy of examining "evolution in reverse" with the four- 
helix bundle ETC hemoprotein cytochrome b562. Start- 
ing with the evolved protein, variants with replacements 
at amino acids near the active site heme were created 
and examined for redox function. We found two general 
results. First, within this conserved protein architecture, 
a range of variation in redox potential e" of about 160 



mV could be obtained within two rounds of (reverse) evo- 
lution, involving only four residues. Statistical analysis 
based on Chebyshev's theorem indicates that this range 
represents, with > 75% confidence, the total range ac- 
cessible through mutations at these positions. Second, 
the wild-type redox potential was not found to be at the 
middle of the chemically accessible range of reduction 
potentials [H, [l3| . Instead, wt b562 exhibits a redox po- 
tential {e° = 167 mV) at the extreme of the chemically 
accessible range (Fig 1). More generally, artificial muta- 
tions on a variety of electron transport proteins of various 
folds and modes of chemical ligation induce redox poten- 
tial changes that span ranges between 100-200 mV (Fig 
1), typically around 150 mV 15, 16, ITJ. Moreover, it is 
possible to sample the majority of the chemically acces- 
sible range through a small number of mutations in the 
vicinity of the active site, with only minimal concomitant 
changes to the structure of the scaffold ^4|. Most im- 
portantly, in nearly every case, these artificial mutations 
push the redox potential in one direction from the wild- 
type value (Fig. 1), indicating that this value represents 
an extremum. In proteins where a few mutations push 
the potential in the opposite direction (e.g., Az. Vin. 
Ferredoxin and Rubredoxin) it is nonetheless clear that 
mutation-induced potential changes are strongly biased 
statistically in one direction from the wild-type potential. 
Maximum likelihood estimation (MLE) of the underlying 
redox potential distributions quantifies this conjecture. 
For instance, in the case of Cyt b562, the nonparamet- 
ric likelihood (see below) that the distribution of redox 
potentials is unbiased is less than 10~*% (Fig. 1). 

III. STATISTICAL CHARACTERIZATION OF 
MUTAGENIC DATA 

Redox potential extremization can be quantified by 
computing the probability of observing a mutation that 
shifts the potential to the opposite side of the putative 
extremum (wild-type, e^t — 0), under suitable paramet- 
ric or nonparametric model distributions. For paramet- 
ric models, the likelihood that the underlying probability 
distribution of the redox potentials obeyed the model dis- 
tributions was assessed by first applyin g th e principles of 
maximal likelihood estimation (MLE) [l^ to determine 
optimal parameters for the model distributions. The em- 
pirical log-likelihood function 

n n 

Z(0;e) = ^log/(£, |0)=^/(0;£,), 

1=1 1=1 

where 9 is the vector of parameters and e is the vector of 
redox potential shifts corresponding to n mutations, was 
maximized by setting the score function vector, defined 
as 

s{e; e) = i{e; = E /(^^ I ^) = E ^(^^ 

4=1 i = l 
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FIG. 1: The shift in redox potential from the wild type value {ewT ~ 0) for active site mutants of several different cytochromes 
and iron-sulfur cluster proteins. Maximum likelihood estimation was employed to quantify the extent to which the proteins 
have evolved toward a redox potential extremum. The likelihood that the true redox potential distribution is unbiased is listed 
below each protein. 



to zero. 

The parametric distributions tested included the nor- 
mal and chi square probability density functions. In the 
case of the normal distribution, analytical expressions 
exist for the MLE estimates; the optimal parameter es- 
timates are equal to the sample mean and sample vari- 
ance. For the normal distribution, the student T-test 
was applied to compute confidence intervals for the mean 
/Lt, and the chi squared test was used to compute confi- 
dence intervals for the variance. For most other distri- 
butions, however, MLE estimates must be determined 
numerically. Two different optimization algorithms were 
used for determining MLE estimates in these cases: 1) 
Newton-Raphson (which employs second derivative infor- 
mation and is the standard MLE optimization algorithm) 
and 2) steepest descent / conjugate gradient. The non- 
central chi square distribution probability density func- 
tion is 

(--1 

fie; ^; so) = Y^i^ - ^or^'^' exp(-(e - eo)/2). 
The parameters v and eg were optimized for the chi 



square distribution. Several initial guesses for these pa- 
rameters were used to seed the MLE optimization algo- 
rithm for the chi square distribution. 

For the optimal normal and chi square distributions, 
the probability of observing a mutation that shifts the 
redox potential to the under-represented side of the wild- 
type potential ( > for Cyt b562, Cyt c, Cyt b5, and 
Rieske ISPs; < for Az. Vin. Ferrodoxin and Rubre- 
doxin) was computed by integration of the probability 
density in this range. The cumulative distribution func- 
tion for the chi square density is 



F(e;i/;£o) = 



7W2,(£-£o)/2) 
r(i./2) 



Tables I and II display the mean, variance, optimal 
MLE parameter estimates, and probability density below 
wild- type for the normal and chi square models. The op- 
timal likelihoods for the normal distributions are gener- 
ally higher than those for the chi square model. However, 
the differences in the chi square versus normal likelihoods 
are significantly lower for the cytochrome proteins. Since 
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the number of parameters in these models is small, and 
because of the skewness of the datasets, it is possible that 
an asymmetric distribution with properties similar to the 
chi square would fit the data better than the normal in 
these cases. 



Figure 2 displays the optimal parametric (normal 
and chi-square) distributions computed by MLE for cy- 
tochromes &562, c, and &5, respectively. As can be seen, 
cytochrome &562 7 which was most extensively mutage- 
nized (and structurally characterized, below) displays 
negligible probability density below the wild type (zero) 
redox potential for the chi square model. For the normal 
model, the integrated density below zero was computed 
for the limiting distributions shifted across the range of 
means corresponding to a 5% confidence interval. Even 
when the variance is upshifted, the probability of observ- 
ing a mutation with redox potential below the mean is 
10.7%, whereas in the opposite scenario this probability 
is only 2.2%, providing a clear indication of redox poten- 
tial extremization in the natural protein. 
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FIG. 2; Optimal parametric distributions computed by max- 
imal likelihood estimation for redox potential shifts. Shifted 
normal distributions correspond to upper and lower limits of 
confidence intervals for a determined by the chi square vari- 
ance test. 

The disadvantage of parametric distributions is that 
the shape of the distributions is constrained by the pa- 
rameterization. Multinomial distributions were used as 
canonical nonparametric distributions. Several nonpara- 
metric multinomial distributions were tested. In this 



case, probabilities were assigned to discrete redox poten- 
tial intervals surrounding the wild-type potential. The 
probability of a mutation falling outside the accessible 
range of 250 mV was set to zero. The interval width 
was set to 250/m or 300/m mV, where m is the number 
of intervals, m was varied between 2 and 8. 

The distributions were scored according to the multi- 
nomial distribution probability density function: 

f{xi, ...,Xm;n,pi, ...,Pra) = — '—yPl^ ' ' ' pI 



where Xi denote the number of observations in interval i. 
Pi denote the respective probabilities, and = "-j 

the total number of observations. 

The likelihoods of several multinomial distributions re- 
sulting from the optimization procedure (for = 6) are 
displayed in Figure 3, for selected proteins. As can be 
seen, even for these distributions, where the probability 
of observing a redox potential in the "forbidden" region 
is generally below 25%, the likelihoods of the models are 
very low - generally an order magnitude below that of the 
optimal model, indicative of extremization of the wild- 
type potential. 



IV. COEVOLUTIONARY DYNAMICS OF 
REDOX POTENTIAL EVOLUTION 

The striking observation of redox potential extremiza- 
tion, confirmed by MLE, begs an evolutionary explana- 
tion. It is clear that there is significant evolutionary se- 
lection pressure acting on the redox potentials; otherwise, 
according to the neutral theory of evolution [l^, they 
would have evolved to maximize robustness to active site 
mutations. Wild-type potentials would then lie in the 
center of the accessible range, and individual mutations 
would alter the potentials by only a small fraction of this 
range - which is not the case. However, there is no ob- 
vious evolutionary advantage to the redox potentials be- 
ing extremized by this selection pressure, since maximal 
fitness (ATP production) follows from maximization of 
the proton concentration gradient, which does not bear 
a simple physical relationship to the redox potentials. 

Direct evolutionary selection for extremized redox po- 
tentials is implausible statistically as well as biophysi- 
cally, based on additional data regarding the distribu- 
tion of potentials within the cytochrome c' family, whose 
members may be perceived as points along a single dy- 
namical evolutionary trajectory. Two out of four mem- 
bers (Cyt c' Chr. Vinosum and Cyt c' R. rubrum) have 
redox potentials at the lower extreme (-5 and -8 mV, re- 
spectively) and two members (Cyt c' Ale. Denitrificans 
and Cyt c' Rps. Palustris) have potentials -1-100 mV 
and -1-130 mV [l^. Hence, it appears that the redox po- 
tentials of naturally-occurring cytochromes are not only 
extremized, but may be alternately maximized and min- 
imized during the course of evolution through a process 
that requires relatively few mutations. Even if evolu- 
tionary selection acted directly on the redox potentials. 
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TABLE I: Statistical properties of mutation-induced redox potential variations in electron transport proteins: 
MLE estimation of normal distribution. "ML" = maximal likelihood. "Conf denotes +/~ confidence intervals for 
mean(/i), standard deviation (a) calculated by the T-test or chi squared test. "Prob forbid" denotes cumulative probability 
density in the "forbidden" region according to the model (e < for the cytochromes and Rieske ISPs, e > for Ferrodoxin 
and Rubredoxin). "Opt" denotes optimal normal distribution; "10%cr" denotes the distributions at the limits of the standard 
deviation confidence interval. For example, for Cyt b562 mutants, the 10% confidence interval for the mean redox potential of 
-60 mV is ±12 mV, and the 10% confidence interval for the standard deviation of 37 mV is ±7 mV. Under the optimal normal 
model, the probability of observing a mutation with redox potential below e = is 5.3%; under the normal model with a — 30, 
this probability is 2.2%, whereas under the normal model with a — 44, this probability is 10.7%. The ML of -186.3 indicates 
that the true redox potential distribution is less likely to be normal for Cyt b562 than it is for the other proteins. Typical error 
bars on redox potential measurements were ±2 — 3 mV. 



^ 0.5 

n 

o 







Likelihood 
1.00 




_Bj ^ 0.5 

.D 





-179-137 


-95 -53 -11 


31 




Likelihood 
1.00 












-175 -125 -75 


-25 


25 75 



C ^ 0.5 

I — ' n 








Likelihood 






1.00 









-79 -37 5 47 89 131 



0.5 



Likelihood 
0.018 




-179-137 -95 -53 -11 31 

Ftedox interval midpt (mV) 



0.5 



Likelihood 
0.233 




-175 -125 -75 -25 25 75 

Ffedox interval midpt (mV) 



0.5 





Likelihood 
0.37 









-79 -37 5 47 89 131 

Ftedox interval midpt (mV) 



0.5 



Likelihood 
0.002 



-179-137 -95 -53 -11 31 



0.5 



Likelihood 
0.108 




-175 -125 -75 -25 25 75 



0.5 





Likelihood 






0.18 










-79 -37 


5 47 


89 131 



FIG. 3: Nonparametric multinomial redox potential distributions for A) Cyt 6562, B) Cyt 65, and C) Az. Vin. Ferrodoxin; 
N=6, with associated relative likelihoods (optimal = 1). Columns 2 and 3 correspond to "2nd, 3rd" distributions from Table 
III. Note the rapid falloff in likelihood with increasing probability density in the "forbidden" region ( > for A,B; < for C; 
redox potentials listed are shifts with respect to wild-type potential). 
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TABLE IL Statistical properties of mutation-induced 
redox potential variations in electron transport pro- 
teins: MLE estimation of chi square distribution. 

"Prob < 0" denotes cumulative probability density in the 
"forbidden" region according to the model. 



it would be necessary to assume the selection pressure 
oscillates due to environmental dynamics that have no 
relation to the known function of the ETC. Such a model 
is not robust to functional form misspecification of the 
fitness measure, and must be rejected if a simpler fitness 
measure requiring fewer extrinsic parameters can explain 
the extremization. 

Optimal control theory provides an explanation for 
the observed behavior that is fully consistent with cur- 
rent evolutionary theory, based on minimal additional 
assumptions. In order to apply OC, it is first nec- 
essary to formulate the dynamical equations governing 
the evolution of the ETC. The terminal oxidation stage 
of the electron transport chain consists of a linked set 
of protein-catalyzed substrate oxidation steps, several 
of which are coupled to protein-catalyzed proton pump 
steps. As described above, electron transfer to the re- 
dox centers alters the pK^j of amino acids involved in 
proton transport and hence indirectly impacts the ef- 
ficiency of the proton pumps. The z-th enzyme acts 
on its substrate through a redox process specified by 
the potential ei(i) as a function of evolutionary time 
t. The role of this i-th enzyme in the fitness measure 
can be characterized by its current evolutionary state 
Xi (i.e., the proton gradient produced by its associated 
proton pump) prescribing the functional utility of the 
enzyme for the energy transduction process. Since the 
efficiency of the proton pumps is a function of the re- 
dox potentials, it is natural to view the network as an 
input-output control system, with the controls consisting 
of e(t) — (ei(i), £2(^)1 ■• • j^Nit)) and the output being 
the system state vector x{t) — {xi{t),X2{t), ■ ■ ■ ,XN{t))- 
Evolution is assumed to be maximizing a biologically 
beneficial function <i>(x) of the chain's state (i.e., the to- 
tal amount of ATP produced) both directly with respect 
to the state x as well as indirectly through the controls 
s{t). 

This evolution of the chain can be modeled in terms 
of the coevolutionary dynamics ^3] of coupled quasis- 
pecies sequence families A and B, corresponding to each 
protein's state and control sequences, respectively. These 
families are described by the multinomial probability dis- 



tributions 

Pa = {uk \ l < k <n^ k"} 
Pb = I 1 < A: < m = k^} 

where k is the monomer alphabet length and /i are 
the respective sequence lengths. The probability of pro- 
ducing sequence Ai as an error copy from sequence A^ is 
given by the elements of the mutation matrix , 

'''''''' 

Wki = Wo (^^^^ j ' (1) 

{fc,^Gl,2,---,K^} where w is the fidelity of (base) repli- 
cation and Wo = w'^ and d{l,k) denotes the Hamming 
distance between sequence k and sequence I (number of 
monomer positions in which they differ). 

The quasispecies kinetic model assumes sequence 
growth by first-order autocatalysis and death by first- 
order decay. We denote by Rk the first-order rate 
constant/parameter for autocatalytic amplification, (i.e. 
replication catalyzed by template Ak, of which fraction 
Wkk leads to identical replica) and by the rate con- 
stant for decay of sequence k. The quasispecies dynam- 
ical equation for the evolution of the probability distri- 
bution PA{t) is then given by 

cik ^^WkiRiai - Dkttk. (2) 

In the quasispecies model, the fitness measure $ enters 
implicitly into the evolution equation through its impact 
on the growth and death rate constants R and D. 

In the coevolution of sequences A and B, the growth 
(and death) rates of the DNA sequence encoding the en- 
tire protein that includes subsequences Ak and Bk' are 
explicit functions of only the state subsequence Ak ■ The 
probability of a control sequence Bk' being replicated 
(or decaying) is then determined by the Ak to which it is 
physically coupled. The probability of subsequences Ak 
and Bk' appearing in the same strand is ak ■ bk' ■ In elec- 
tron transport proteins containing both a proton pump 
and a redox center, these correspond to the pump and 
active site sequences. We denote the mutation matrix 
for the second sequence by and the growth and death 
rate constants by S and E. Then, in the quasispecies 
evolution equation for subsequence B, the growth rate 
constant is completely determined by the matrix W, the 
vector R, and the probabilities ak- 
in coevolutionary quasispecies dynamics, the effect of 
a second coevolving species on the fitness measure $ of 
the first is typically modeled through a perturbation of 
the first-order rate constants R, D, or both [23|- There- 
fore, in accordance with subsequence B functioning as 
a control, assume that sequence B can perturb the fit- 
ness measure $ such that it affects the growth rate con- 
stant of A; i.e., introduce a 6-dependence in Ri, writing 
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Ri{bi, ...,bm) or Ri{Pb)- We then have: 

n 

flfc = '^WMRi{PB)ai - Dkak, 



1=1 



Y^VuSih 



kOk 



1=1 



1=1 



kl 



.q=l \p=l J 



2^ Dpflp I hk- 
\v=i ) 



We assume the effect is restricted to the growth constant 
R. It is natural to work within a first-order model where 
Riibi, bm) is a linear function of the fo^'s. In the first- 
order approximation, Ri{bi, ...,bm) can be written i?° + 
J2^=i R'lr^r, such that we get 



1=1 



Associated with each state sequence Ak is the value 
-Ffe G R of a component of the associated physical state 
vector X of the protein network (respectively Hk for the 
control vector e). The expected values of the compo- 
nents of the state and control vectors of the protein net- 
work are then Xi = {xi) = J2k=i ^k^ o-k^ ^ = i^i) ^ 
J2T=i ^k^'^^k ■ '^^^ tertiary structure of the protein mi- 
croenvironment surrounding the redox center con- 
strains Ei (t) to a finite range 



(3) 



Because protein tertiary structure is less flexible than sec- 
ondary structure, it is reasonable to assume the bounds 
s\{t) and sf{t) vary at a slower rate than the redox po- 
tential £i{t) during evolution. 

The evolution of the expectation value of the state 
corresponding to evolution of the distribution PA{t) = 
{ak{t) I 1 < fc < n} is then given by 



djx^jt)) 
dt 



k 



k=l 



fc=i 



E ^kiRkibi, ■■■ , bm)ak - D 



kO-k 



.1=1 



n n 



= E ^'^ E ^kiRiibi, • • • , b^)ai - J2 DkFkak 

_l 



k=l 1=1 



k=l 



Now, because the multinomial distribution Pg is sharply 
peaked with a small variance around a master sequence 
"^max (see below), such that (e^ 
sonable to make the replacement 

R° + -R>max = R° + R'falH,,_____ 

copy fidelity w « 1, the off-diagonal elements Wij{ i ^ 
j) Wii. Under these approximations, we can write: 



-^max^max; it is rca 
J ' J2r=lRjk^k ~ 

Furthermore, since the 



d{x,{t)) 



n n 



dt ' Y^kY^kiiR'i +R'i£a/H^,^)ai 

k=l 1=1 

n 

- DkFkak- 



k=l 



If interactions are permitted between state vector com- 
ponents Xi{t), this can be written compactly as: 



dxi{t) 
dt 



/,(x(i),i)-|-g,(x(t),t)£,(t). (4) 



V. OPTIMAL CONTROL OF EVOLUTIONARY 
DYNAMICS 

We note that the above evolutionary dynamics frame- 
work is based solely on the quasispecies theory and the 
known function of the ETC. We now show that the ob- 
served extremization implies that the rate parameters 
R'l^ and DJ,^ have been set such that the e(i)'s are op- 
timal for maximizing the increase in evolutionary fitness 
in a given evolutionary time step dt. This entails a maxi- 
mization of $(x) with respect to the controls e{t), subject 
to the inequality constraint in Eq. ([3]) and the dynami- 
cal constraint in Eq. It is convenient to rewrite the 
inequality constraint in the form of an equality through 
the introduction of so-called slack variables [t) where 



Gi{t) — Gi{ei, 



6 ■ S- 



(5) 



\{t)){eUt)~e^{t))-edt)^^■ (6) 



When each of these slack variables [t) is allowed to take 

on arbitrary real values, then the equality constraint in 
Eq. ([5]) is consistent with Eq. ([3]). We may now define 
the fitness measure J as having the following form: 



j=<i>(x)+^ r mG^{t)+Y r^^w 



■^Xi - ft - 9i£i{t) 



dt. 



(7) 



The introduction of the Lagrange multiplier functions satisfied, respectively. Equation ([7]) leads to the biolog- 
\i{t) and (ii{t) will assure that Eqs. (jl]) and ([5|) are 
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FIG. 4: (a) The evolution of the redox potential ei{t) for the i-th enzyme within optimal control theory. The early evolutionary 
period near t « is unspecified. During evolution, the potential "bangs" from its lower and upper accessible values, e'(t) and 
ef{t), respectively, at critical times ti,t2, • ■ ■ where effective mutations have occurred. The current evolutionary time is T. (b) 
The evolutionary time dependence of the product Xigi of the Lagrange function Ai and the control coupling function gi. The 
zero crossings of Xigi occurs at ti,t2, ■ ■ ■ where the redox potential ei{t) undergoes evolutionary jumps in (a). 



leal evolutionary process expressed as max^(t) J. Max- 
imization of J can be treated as a problem in the cal- 
culus of variations, with the unknown functions being 
the elements of the vectors e, /3, ^, x, A. A variation of J 
with respect to these functions will produce a set of non- 
linear equations whose solution would specify the state 
of the evolving protein network from its initial condition 
at t = to the current time T. Since we have not com- 



pletely specified the functions /i(x(t), t) and gi{x{t),t) in 
Eq. a detailed study of the evolutionary dynamics 
cannot be carried out here. However, for our purpose 
of analyzing the mutation data above, we do not need 
this level of detail. It is sufficient to only consider varia- 
tions of J with respect to ^, /3, and e, which produce the 
following equations 



5J _ 
5J 



-2A(i)6(t) =0 
= G^{t) = 0. 



5J 



^(3,{t) [-2e,{t)+e\[t)+e't[t)\+\[t)g,{t)=Q 



(8) 
(9) 
(10) 



We may now analyze the evolutionary consequences ^^(t) or Pi{t) = 0. Considering the first case, £,i{t) = 
of these equations. First, Eq. ([5]) implies that either 



9 



0, it is evident from Eqs. ^ and (O that the redox 
potential ei{t) must take on the value ei{t) = e"(t) or 
El (t) = e[ (t) . We may then solve for /3j (t) from Eq. (fTI 
by first defining di as 



d,: = - 



such that 



2e,{t) + slit) + snt) (11) 



(3,{t) = -K{t)g,{t)/d,{t). (13) 



The second circumstance, Pi(t) — 0, implies that ^i{t) is 
free to take on any value prescribed by Eq. (O, given 
that ei(t) is restricted to the domain specified in Eq. 
([3]). In this case, it is also evident from Eq. (fTO|l that 
Xi{t)gi{t) — 0, which is expected to only be valid at dis- 
crete times t — tn, n — 1,2, ■ ■ ■ . These time points t„ 
denote the locations where the control field "bangs" from 
one extreme limit of the range to the other in Eq. ^ 
during evolution. 

This behavior may be explicitly seen by considering 
the curvature 



summit'] 



= -2(3,{t)5{t - t') < 



(14) 



where 5{t — t') is a Dirac delta function, and the inequal- 
ity corresponds to requiring that J be maximized. Thus, 
for the case £i{t) = £"(^) in Eq. (fTTj) . it follows that 
di < 0, thereby corresponding to Xi{t)gi{t) > 0, to as- 
sure that Eq. (fTi]) is satisfied. Similarly, in the oppo- 
site case of ei{t) — e\{t), we have that di > and that 



Xi(t)gi{t) < 0. The points t„, n = 1,2, ■ ■ ■ correspond 
to the times at which Xi{t)gi(t) changes sign by passing 
through zero. This behavior is indicated in Figure 2. The 
possible evolution of the extremum values {t) and e" [t) 
is also indicated in the figure. 

Importantly, the redox potential data above [l3l. O[l5l 
El; [13 ^'^'^ fully consistent with this analysis of bang-bang 
control behavior. That is, at the present evolutionary 
time T, each redox potential Si (T) should be at a locally 
accessible extreme value. The introduction of artificial 
mutations in the laboratory could then only take £i{T) 
away from its extreme value in a consistent direction for 
each protein, as found above. Moreover, assuming mem- 
bers of the cytochrome c' family lie along the same evo- 
lutionary trajectory, their alternatively maximized and 
minimized redox potentials are consistent with the above 
model for t < T. We emphasize that this finding of 
optimality is based solely on statistical inference and 
variational calculus and does not imply anything about 
the mechanism by which optimality is achieved. How- 
ever, the required tuning of the rate constants R'lr-i^'ks 
to optimal constant values is straightforward to achieve 
via reorganization of the protein's tertiary structure 14 1 
through genetic recombination, and avoids the biophysi- 
cally implausible assumption of direct evolutionary selec- 
tion for redox potential extremization on an oscillating 
fitness landscape. 

The OC prediction of bang-bang control behavior is 
contingent upon the circumstance that the cost func- 
tional does not explicitly depend on the controls (except 
through the Lagrange multiplier that imposes the dy- 
namical constraint), such that 



J 



J = $(x) + E A(t)Gi(t) + E ^i(t) [^^' - - 



(15) 



r 



If auxiliary penalty terms C explicitly depending on the 
controls ei{t) are introduced, then J ^ J — C, and the 
optimal controls need not be singular, i.e., they may not 
abruptly "bang" from one extreme to the other. 
Two biologically plausible scenarios correspond to: 



C 



'E 



, d£i{t) s 
' dt ' 



dt. 



(16) 



which places a cost on the rate at which control changes 
occur, and 



. Jo 



(17) 



which places a penalty on the time-average of the control 
magnitude, in addition to restricting this magnitude to 



a bounded range. 

Figure 3 displays possible optimal controls £{t) (redox 
potentials for the ETC) resulting from these respective 
cost functionals. In the case of the former cost, bang- 
bang control can still be produced, but with a rounding- 
off of the sharp corners at the jump times. Biologically, 
this corresponds to only having an ability to make evo- 
lutionary changes in a gradual fashion, while still tak- 
ing advantage of the extreme accessible controls as being 
biologically most effective. By contrast, the latter cost 
may allow the controls to take on any intermediate val- 
ues over evolutionary time, since extreme control mag- 
nitudes are highly penalized. The former scenario has 
a natural structural interpretation, since a given change 
in a functional protein property, such enzyme activity, 
may require multiple (active site) mutations occurring in 
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FIG. 5: (a) The evolution of the redox potential ei{t) for the i-th enzyme in the case that the evolutionary cost functional is 
J' — J — C, where J is given by Eq. (7) and C by Eq. (8). The added influence of the cost in Eq. (2) builds in the ability 
to make smooth evolutionary changes (red curve) in the redox potential while still remaining within its accessible lower and 
upper values of £i{t) and e^{t). The dashed line corresponds to just operating with the cost functional in Eq. (7) (see Figure 
2(a) FIX), (b) The cost functional is J' = J — C, with C given by Eq. (9). In this case, the genetic pressure on the magnitude 
of the redox potential permits it to smoothly evolve to have any value within the accessible range between £\{t) and e"(i). 



succession rather than in unison. The quasispecies error 
threshold sets a limit on the number of mutations that 
can be borne by an evolving population per generation 
[2l| . Although bang-bang behavior may not be as appar- 
ent in such cases, optimal control may still be in effect. 

VI. CONCLUSION 

A natural question concerns the generahty of optimal 
control phenomena in evolutionary dynamics. Optimal 
control could in principle be operational in any system 
where evolution of the central function of a protein net- 
work is coupled to the evolution of an ancillary protein 
function. Our results indicate that it is worthwhile to 
investigate whether the evolutionary dynamics of other 
biochemical networks with coupled functions exhibit the 
characteristic signatures of being under optimal control. 
Bang-bang extremization, while not the only such sig- 



nature, is simple to detect and provides compelling evi- 
dence for underlying OC phenomena. Such optimal con- 
trol strategies have a particularly natural interpretation 
within the general framework of evolutionary optimiza- 
tion. 

The observation that coevolving biopolymer sequences 
may optimally control each other's evolution raises the 
prospect of artificial optimal control of evolutionary dy- 
namics. Possible applications include the control of repli- 
cation fidelity in nucleic acid amplification reactions and 
the design of therapeutics that dynamically regulate the 
evolution of viral populations. 
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